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Abstract 

We carry out a microscopic analysis of the ground and excited states of the Nas 
metal cluster within the jellium model. We perform a series of configuration inter- 
action calculations on a Hartree-Fock basis and construct eigenstates of the Hamil- 
tonian which carry up to 4-particle 4-hole components. Based on the analysis of the 
dipole transition strengths, we single out those states which can be interpreted as 
the collective dipole plasmon and its double excitations. These modes are found to 
possess a high degree of harmonicity, deviations from the harmonic limit remaining, 
however, of the order of 10%. 
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Collective vibrational states are known both in metal clusters [1,2] and in nu- 
clei [3]. They are interpreted as the excitations of vibrational quanta, the so 
called plasmons or phonons. In particular, the dipole plasmon excitation in 
clusters, corresponding to an oscillation of the centre of mass of the valence 
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electrons against that of the positive ions, is a very collective mode and dom- 
inates the response to a laser field. Microscopically, its frequency and strong 
collectivity are quite well reproduced by the random phase approximation 
(RPA) within the jellium model. This model [4,5] considers a uniform positive 
charge distribution generated by the ions (the jellium) which interacts with a 
cloud of delocalized valence electrons (also interacting among themselves) via 
the Coulomb interaction. RPA predicts that the plasmon oscillation is per- 
fectly harmonic [6], i.e. that states corresponding to a n-fold excitation of the 
plasmon exist and their energy is equal to n times that of the single plasmon. 
If the jellium-electron interaction can be approximated by a harmonic oscil- 
lator potential, as it is reasonably true for a highly positively ionized cluster, 
this is what one expects on general grounds [7] . Indeed, in that case the total 
electron eigenfunctions are simply products of the type 

(ri, r 2 , ...r N ) = ^(R)^^, r' 2 , ...r' N ) , (1) 

where R is the coordinate of the centre of mass of the electrons with the har- 
monic oscillator quantum numbers (nl) describing its motion, while r[ are the 
intrinsic coordinates of the electrons with v and A classifying the harmonic 
oscillator states of their relative motion. A dipole external field acting on the 
ground state of the system, t/>oo0oo, causes the transition to ^oi^oo, acting 
on the latter excites the states t/>io0oo an d Vte^oo, and so on, generating a 
perfectly harmonic band. Of course, other bands exist, based on different in- 
trinsic motion states 0„a, with no electromagnetic transitions among them. 
This was checked numerically in Ref. [8] for the very simple case of two inter- 
acting electrons moving in a harmonic oscillator potential. When this potential 
is replaced by the Coulomb potential generated by the jellium charge distri- 
bution the above level scheme is modified and it was found in Ref. [8] that 
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different bands still exist, deviating from the harmonic limit and with non 
zero multipole transition probabilities among them. These results show that 
the anharmonicities in the dipole plasmon excitation are due to the coupling 
between the intrinsic and centre of mass motions of the electrons. 

This has been recently discussed quite in general in Ref. [9]. In this work, the 
authors make explicit in a very clear way the coupling between intrinsic and 
centre of mass motions by considering an expansion of the ionic background 
potential V ion (ri) = l / ion (r l '+i?) in power series with respect to R. Appreciable 
deviations from the harmonic scheme are found within this approach [9]. For 
example, for a Nag 3 cluster it is found that the most important optical dipole 
transitions are separated by a AE n equal to 2.8 eV, 3.07 eV and 3.26 eV for 
the single (n = 1), double (n = 2) and triple (n = 3) plasmon excitations, 
respectively. These results should be compared with the contradictory ones 
found in Refs. [10] and [11], namely essentially zero and huge anharmonicities, 
respectively. This is due to the fact that the coupling between intrinsic and 
centre of mass motions is almost completely neglected in Ref. [10]. On the other 
hand, the boson expansion used in Ref. [11] was truncated at the same order 
as previously done for the study of the double excitation of Giant Resonances 
in atomic nuclei [12]. In the latter case, anharmonicities of the order of a few 
hundred keV as compared with the harmonic limit energy of 20-=-30MeK were 
found. The huge anharmonicities found in metallic clusters are probably an 
indication that the convergence of the boson expansion is much slower in that 
case and this can be related to the long range of the Coulomb interaction [7]. 
To our knowledge, there is no experimental clear evidence on the existence of 
states corresponding to the double excitation of the dipole plasmon [13,14]. 
On the other hand, in atomic nuclei, the existence of multiphonon states has 
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been known for many years [3] and the anharmonicities in their excitation 
spectra have been found to play a role in several physical processes [15]. 

This and the very interesting results found in Ref. [9] encouraged us to attempt 
a configuration interaction (CI) calculation for small metal clusters. Indeed, 
this is a priori the best method to obtain the energy spectrum. On the other 
hand, since the states we want to study are quite high in energy, namely at 
about twice or more that of the plasmon, the configuration space required 
to get reliable solutions becomes very rapidly prohibitively large when the 
size of the cluster increases. We have thus limited our analysis to Na 8 . Similar 
calculations for several light Na clusters have been reported in Ref. [16] where, 
however, only the ground state and the singly excited dipole state were studied. 

Let us then shortly describe our calculations and discuss the results. Within 
the jellium model the motion of the valence electrons is determined by the 
hamiltonian 



with 



and 



# = + (2) 

i i<j 



(l/2r c )(r 2 /r c 2 -3) forr<r c 

Ze 2 { 



— 1/r for r >r c 

where r c is the radius of the jellium sphere, i.e. r c = r s N^ 3 with r s the Wigner- 
Seitz radius which is 4 a.u. for Na and N e is the number of electrons. We are 
aware of the fact that, for a small metal cluster, the jellium approximation may 
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be not completely adequate for a quantitative comparison with experimental 
data. However, it already contains many important physical features while 
allowing not too heavy calculations as compared with more elaborated models. 

As a first step we make a Hartree-Fock (HF) calculation in order to fix single 
particle (sp) energies and wavef unctions. This is done by allowing the wave- 
functions to be superpositions of harmonic oscillator wavefunctions (h.o.w.f.'s). 
The number of h.o.w.f.'s has to be chosen large enough to get a HF ground 
state energy satisfactorily stable for small variations of the h.o. parameter 
around the value giving the minimum. Next, we construct all Slater determi- 
nants with fixed values of the projection Ml of the total angular momentum, 
of the projection Ms of the total spin and of the parity n. This set is truncated 
by i) truncating the sp basis; ii) putting a maximum value n for the number of 
particle- hole excitations (np-nh) with respect to the HF ground state and/or 
for the unperturbed excitation energy. We have considered several truncations 
corresponding to i) up to 10 HF orbitals above the Fermi level; ii) all Slater 
determinants having up to 3p-3h configurations and containing those Ap-Ah 
configurations whose unperturbed energy with respect to the HF ground state 
is less than a given cutoff energy E c . The use of the HF basis should optimize 
the convergence of the results because part of the correlations are already 
taken into account in the reference state. 

The largest basis we have been able to manage has dimension ~ 700000; it 
is to be noted, however, that many matrix elements of the hamiltonian in 
such a basis are zero: the hamiltonian matrix is sparse. We have then used a 
NAG library routine especially intended for such a case. The routine finds the 
N eigenvalues of largest absolute value and the corresponding eigenvectors. 
This method is very suitable for our case since we are interested in the lowest 
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negative eigenvalues which are the largest ones in absolute value. The time 
required by the routine to find the solutions increases very rapidly with N. 
Therefore, N has to be taken as low as possible. From RPA we know that the 
excitation energy of the dipole plasmon is located around 3 eV. Since we are 
interested in single and double dipole plasmon excitations, we have to choose TV 
such that all eigenvalues up to ~ 7 eV are determined, and this would mean 
a quite large N since the states of the basis we use do not have a definite 
value of the total angular momentum and spin. On the other hand, since the 
hamiltonian commutes with L and S, its eigenvectors are also eigenvectors of 
L 2 and S 2 in addition to L z , S z and parity. We take advantage of this and 
use the following procedure to select the eigenstates of H belonging to some 
definite value of the angular momentum L. By running the calculation with 
M L = L, the eigenstates with L < L are trivially eliminated. By adding to the 
hamiltonian a term a[L 2 — L(L + 1)] with a > the eigenvalues corresponding 
to L > L are shifted up. Therefore, if a is chosen large enough, only those with 
L = L are selected as the lowest ones (i.e. with the largest absolute values). 
In reality one must be careful since, if a is chosen very large, the eigenvalues 
corresponding to L ^> L are so much shifted up that they become positive 
and the largest ones. Therefore, a has to be kept not too large and then 
a few eigenvalues not corresponding to L may be mixed in the region of the 
spectrum we look at. However, this problem is easily eliminated by calculating 
after the diagonalization the angular momentum of each eigenstate. The same 
reasoning applies also for the spin and we add to the hamiltonian another term 
(3[S 2 — S{S + 1)] with (3 > 0. The above sketched procedure is very effective 
and it has allowed us to limit the number of eigenvalues we are interested in 
within a maximum value of 50, for each of the three cases relevant for our 
scopes, namely U S = 1~0, + and 2 + 0, and an excitation energy less than 
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7 eV. The energies of the lowest 12 HF sp states are reported in Table 1. 

Each sp state has been expressed as superposition of h.o.w.f.'s, with principal 
quantum number running from to 8. The HF states are labeled by the 
principal quantum number n of the predominant component and the angular 
quantum number /. 

For a start, we have performed two series of calculations by considering those 
sp states and including all Slater determinants with up to 2p-2h and 3p- 
3h configurations. By comparing the two series of results we have observed 
that the 2p-2h space is far from being sufficient for a good description of the 
states we are interested in. Indeed, the inclusion of the 3p-3h configurations 
strongly modifies the energies of the excited states and, to a less extent, of the 
ground state (see Table II). In order to go further one has to include Ap-Ah 
configurations. However, in this case the number of configurations is too large 
and it is necessary to introduce an energy cutoff E c , i.e. to include only those 
Ap-Ah Slater determinants whose unperturbed energies are not higher than E c 
above the HF ground state. By repeating the calculations with E c increasing 
from E c = 19 eV to E c = 25 eV we get a lower and lower ground state energy 
as shown in Table 2. Looking at the third column of the table, where we show 
the difference between the ground state energies relative to two successive 
calculations, one can conclude that a very good numerical convergence has 
been reached for the ground state energy. Indeed, the values obtained with 
E c = 24 and 25 eV differ only by 2 meV. However, we are also interested in 
excited states having a quite high energy and convergence must be checked 
also for them. In order to do that we have calculated the root mean square 
value a 

= /WEE A? of the shifts A; obtained in two calculations with E c 
differing by 1 eV. By considering all states having excitation energy less than 
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7 eV, for the case L n S = + 0, we have got a = 0.064 eV when going from 
E c = 20 eV to E c = 21 eV while tr = 0.039 eV in going from E c = 24 to 
E c = 25 eV. It is also to be noted that, in the latter case, the largest shift is 
A = 0.046 eV corresponding to 0.7 %. Similar results have been obtained for 
the 2 + S = and l~S = spectra. 

A further comment has to be added. In order to be able to increase so much 
the energy cutoff for 4p-4h configurations we have followed a suggestion of Ref. 
[16]. Namely, the high angular momentum Og and Oh single particle states have 
been suppressed from the basis. We have checked that, for E c = 20 eV, the cal- 
culations with and without those states give almost indistinguishable results, 
more precisely a maximum shift of 0.2 % in the considered energy region. The 
next two sp states above the twelve ones we have considered in the basis have 
1=6 and 1=7. Therefore, their inclusion should not modify the spectra. Still 
above there is a 2d orbital, but its HF energy is 2.19 eV and we can reason- 
ably neglect its contribution. From the above considerations we conclude that 
a satisfactory convergence has been reached and the so obtained energies and 
wavefunctions can be reliably used in the analysis of anharmonicities. 

In order to single out the states which can be interpreted as the collective 
dipole plasmon and its double excitations we have calculated the electric dipole 
transition strengths from the initial states \ipi >, + S = and 2 + S = 0, 
to the final states \ipf >, 1~ S — 0. More precisely, we have calculated the 
following quantity: 

\T fi (El)\ 2 = \^C ] <^fW a a a/ \^>\ 2 , (5) 

aa' 
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with a = (n a , l a ,Ca) and 

= (« II Yl II «') = Rn a l a ,n a ,l a ,(L II ^1 II kO$r«<r a , , (6) 

Rni,n>i> = J r 3 ip* nl (r)ip nn/ (r)dr . (7) 

The (/ || Y\ || /') in eq. (6) are the reduced matrix elements defined as in Ref. 
[15] and ip n i in eq. (7) are the Hartree - Fock sp wave functions. In Fig. (1) 
we report the largest \Tfi(El)\ 2 values (greater than 0.8 A 2 ) obtained in the 
calculation performed with 12 HF states and within the Ap-Ah space with E c 
=25 eV. As we are interested in the study of the collective dipole plasmon 
and its double excitations, we show only the transitions involving the 1~ states 
located around the single plasmon excitation energy. 

The selected levels can be grouped in two "bands" based on the ground state 
and on the lowest l - S = state (If), respectively. The two bands are essen- 
tially not connected by dipole transitions. The second 1~ S — state (lj) at 
3.07 eV excitation energy is strongly coupled to the ground state and can be 
identified with the dipole plasmon excitation. The only three excited states 
having a large dipole transition strength to this (lj) state are two 2 + S = 
states, at 6.41 and 6.42 eV, and one + S = state at 6.60 eV, corresponding 
to an energy jump of 3.34, 3.35 eV and 3.53 eV respectively, to be compared 
with the 3.07 eV excitation energy of the (I2 ) state. They might be identified 
as corresponding to double excitations of the dipole plasmon. Indeed, they 
are very close to each other and their energy is not far from twice that of the 
single plasmon. The deviations from the harmonic limit are of the order of 
10%. As can be seen from the figure, the strength of the transition connecting 
the (I2 ) state to the two-plasmon states is 12.12 A 2 to be compared with the 
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harmonic value of 12.36 A 2 corresponding to the double of the one-plasmon 
strength. We also note that this strength is fragmented between the 2 + S = 
and + S = states in a proportion close to the harmonic limit, i.e. 4/3 and 
2/3 respectively. 

The 2 + S = and + S = having a strong transition to the (If) state 
(see the "lateral band" in the figure) lie at 3.40 eV and 3.65 eV above it. 
This makes plausible their interpretation as collective dipole excitations built 
on top of the (lj~) state, all these states being characterized by an intrinsic 
motion wave function </> 01 (see eq. (1)). This also justifies the extremely small 
value (0.002 A 2 ) of the transition strength between the ground state and the 
(li) state. 

The quality of the results obtained in the present calculations can be judged 
by looking at sum rules. It is well known [6] that, if |0) and \v) are the exact 
ground and excited states of a system, then the following equality holds 

^(E u -E )\(u\T x \0)\ 2 = l(u\[T x ,[H,T x ]]\0) , (8) 

where H is the total hamiltonian and T A the transition operator of multipo- 
larity A. The r.h.s. of the above equation can be evaluated exactly and is (see 
eq. (2.47) of Ref. [15]) 

EWSR = ^ A (2A + l) 2 2A _ 2 

where m is the mass of the electron and (r 2X ~ 2 ) is the expectation value of 
the indicated quantity in the ground state. Since we are looking at A = 1 
transitions, the EWSR is completely independent of the ground state and is 
equal to 21.83 A 2 ■ eV. The l.h.s. of eq. (8) turns out to be 21.62 A 2 ■ eV which 
is 99.04% of EWSR. Therefore, we can conclude that our calculated energies 
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and wavefunctions satisfy up to an extremely good level this very stringent 
condition. In order to further check the numerical convergence we have com- 
pared the results obtained with 12 HF orbitals and cutoff energies equal to 
E c = 22, 23, 24, 25 eV , finding a smaller and smaller variation at each step, 
the last one being 0.06 A 2 ■ eV, i.e. 0.2 %. It is also worth mentioning that 
the contribution of the second \~ state to the sum is 19.01 A 2 ■ eV, which 
means 87% of the total. Therefore, its identification with the collective dipole 
plasmon state is very well justified. 

By summarizing, in this work we have carried out a microscopic analysis of 
the ground and excited states of the Nas metal cluster within the jellium 
model. We have performed a series of configuration interaction calculations 
in a Hartree-Fock basis and constructed eigenstates of the Hamiltonian which 
carry up to 4-particle 4-hole components. Based on the analysis of the dipole 
transition strengths, we have singled out those states which can be interpreted 
as the collective dipole plasmon and its double excitations. These modes have 
been found to possess a high degree of harmonicity, deviations from the har- 
monic limit remaining, however, of the order of 10%. These values are consis- 
tent with the anharmonicities found by Gerchicov et al. [9], using a different 
technique, in heavier clusters. We want to stress that in principle a config- 
uration interaction calculation is the most accurate approach to reproduce 
the spectrum of the system and then to evaluate its anharmonicities. On the 
other side, a limitation of this kind of calculations is that they are very heavy 
numerically and this is the reason why we limited our analysis to the small 
cluster Na 8 . 
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Fig. 1. The largest \T fi (El)\ 2 from the 0+ S = and 2+ S = states to the 1" 
S = ones. The values, in A 2 , are indicated in brackets. The calculation is done 
within the to 4p-4h space with an energy cutoff equal to 25 eV and 12 HF orbitals. 
The states are grouped in two bands based on the ground state and on the lowest 
1~ S = state. For the excited states we report the energies corresponding to the 
transitions indicated by the arrows. The two arrows in correspondence with the 
highest 2 + state refer to two almost degenerate states which are not distinguishable 
in the figure (see the text). 
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0.15 


0.25 


0.50 


0.83 


n I 
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2 


4 
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2 1 


1 3 


E (eV) 


0.86 


1.12 


1.16 


1.46 


1.62 


1.76 



Table 1 

Energies (eV) of the lowest 12 HF sp states; n is the principal quantum number of 
the predominant component and / is the angular quantum number. 
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E (eV) 


A (eV) 


2p-2h 


-142.095 




3p-3h 


-142.176 


-0.081 


4p-4h, E c = 19 eV 


-142.235 


-0.059 


4p-4h, E c = 20 eV 


-142.248 


-0.013 


4p-4h,E c = 21 eV 


-142.252 


-0.004 


4p-4h, E c = 22 eV 


-142.258 


-0.006 


4p-4h, E c = 23 eV 


-142.275 


-0.017 


4p-4h, E c = 24 eV 


-142.283 


-0.008 


4p-4h, E c = 25 eV 


-142.285 


-0.002 



Table 2 

In the second column the ground state energies are reported for each calculation. 
In the third column Ao represents the difference between the corresponding ground 
state energy and that reported in the previous line. 
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